1 Introducción

1.1 Infraestructura informática para el análisis

Este estudio se llevará a cabo usando R/Bioconductor por lo que es preciso tener instaladas un conjunto de librerías. Esto puede hacerse siguiendo el procedimiento descrito a continuación.

El código que se presenta debería ejecutarse tan sólo una vez!

if (!require(BiocManager)) {
    install.packages("BiocManager", dep = TRUE)
}

installifnot <- function(pckgName, BioC = TRUE) {
    if (BioC) {
        if (!require(pckgName, character.only = TRUE)) {
            BiocManager::install(pckgName)
        }
    } else {
        if (!require(pckgName, character.only = TRUE)) {
            install.packages(pckgName, dep = TRUE)
        }
    }
}
installifnot("limma")
installifnot("edgeR")
installifnot("org.Hs.eg.db")
installifnot("clusterProfiler")
installifnot("dplyr", BioC = FALSE)
installifnot("gplots", BioC = FALSE)
installifnot("ggvenn", BioC = FALSE)
installifnot("pheatmap")
installifnot("prettydoc", BioC = FALSE)
installifnot("ggnewscale", BioC = FALSE)

Puede resultar útil, aunque no es para nada imprescindible, disponer de, al menos dos, directorios específicos: - datos (o un nombre similar) donde guardar y de donde cargar los datos. - resultso (o un nombre similar) donde escribir los resultados.

2 Los datos para el análisis

2.1 Lectura de la matriz de contajes

Los datos de contajes se encuentran descargados desde GEO al archivo Rawcounts.csv desde donde se leeran.

Con la información sobre los grupos u otras covariables o variables auxiliares se creará un objeto (el habitual “targets”) que debe estar sincronizado con el anterior.

3 Preprocesado de los datos

3.1 Estandarización de los contajes

Además de filtrar, es bueno expresar los contajes en “CPMs” es decir “counts per million”, lo que no modificará los resultados del filtraje, pero estandarizará los valores, lo que es útil y necesario para los análisis posteriores.

                COV145 COV147 COV149 COV150 COV151 COV152
ENSG00000223972      0      0      0      0      0      0
ENSG00000227232      1      0      3     16      1      2
ENSG00000278267      3      1      2      0      3      0
ENSG00000243485      2      0      0      1      0      0
ENSG00000284332      0      0      0      0      0      0
                    COV145     COV147    COV149     COV150     COV151    COV152
ENSG00000223972 0.00000000 0.00000000 0.0000000 0.00000000 0.00000000 0.0000000
ENSG00000227232 0.06644967 0.00000000 0.2218932 0.84819243 0.07320991 0.1506093
ENSG00000278267 0.19934902 0.07170342 0.1479288 0.00000000 0.21962973 0.0000000
ENSG00000243485 0.13289935 0.00000000 0.0000000 0.05301203 0.00000000 0.0000000
ENSG00000284332 0.00000000 0.00000000 0.0000000 0.00000000 0.00000000 0.0000000

Una vez los datos estan como CPMs, se procede a filtrarlos,

3.2 Filtraje de genes poco expresados

Los genes con recuentos muy bajos en todas las librerías proporcionan poca evidencia de expresión diferencial por lo que es habitual eliminar aquellos genes que, o bien son poco variables, o bien presentan poca o nula expresión en la mayoría de las muestras.

En este caso, siguiendo las indicaciones proporcionadas, se opta por conservar únicamente aquellos genes que presentan algún valor en, al menos, tres muestras de cada grupo.

[1] 60683    34
[1] 33039    34

Aunque no sea más que un ejemplo basta con ver los dos primeros genes para comprobar como el primero no cumple la condición, en el grupo “SANO”, mientras que los siguientes sí que la cumplen. Por lo tanto, al filtrar desaparece el primer gen de la matriz filtrada, pero no los dos siguientes.

                    COV145     COV147    COV149     COV150     COV151    COV152
ENSG00000223972 0.00000000 0.00000000 0.0000000 0.00000000 0.00000000 0.0000000
ENSG00000227232 0.06644967 0.00000000 0.2218932 0.84819243 0.07320991 0.1506093
                  COV153     COV154    COV155    COV156    COV157    COV175
ENSG00000223972 0.000000 0.00000000 0.1163313 0.0000000 0.1173623 0.0000000
ENSG00000227232 1.388085 0.41957493 1.1633130 0.0000000 0.1173623 0.0000000
                   COV176    COV177    COV178    COV179    COV180    HEA061
ENSG00000223972 0.0000000 0.0000000 0.0575402 0.0000000 0.0000000 0.0000000
ENSG00000227232 0.9503052 1.9885554 0.1726206 0.1462825 0.6558384 0.3528731
                   HEA062    HEA063     HEA064    HEA065     HEA066    HEA067
ENSG00000223972 0.0000000 0.0000000 0.06183451 0.0000000 0.00000000 0.0000000
ENSG00000227232 0.0000000 0.0000000 0.06183451 0.1029819 0.05475139 0.0000000
                    HEA068     HEA069     HEA070    HEA071    HEA181    HEA182
ENSG00000223972 0.00000000 0.00000000 0.00000000 0.0000000 0.0000000 0.0000000
ENSG00000227232 0.06713275 0.06693886 0.04887949 0.1122282 0.2003508 0.2800970
                    HEA183   HEA184    HEA185    HEA186
ENSG00000223972 0.00000000 0.000000 0.0000000 0.0000000
ENSG00000227232 0.05465108 0.000000 0.0000000 0.4151761
 [ reached getOption("max.print") -- omitted 4 rows ]
                    COV145     COV147    COV149     COV150     COV151
ENSG00000227232 0.06644967 0.00000000 0.2218932 0.84819243 0.07320991
ENSG00000278267 0.19934902 0.07170342 0.1479288 0.00000000 0.21962973
                    COV152    COV153     COV154   COV155    COV156     COV157
ENSG00000227232 0.15060933 1.3880853 0.41957493 1.163313 0.0000000 0.11736232
ENSG00000278267 0.00000000 0.0000000 0.06992916 0.174497 0.1279917 0.11736232
                   COV175    COV176    COV177    COV178     COV179    COV180
ENSG00000227232 0.0000000 0.9503052 1.9885554 0.1726206 0.14628246 0.6558384
ENSG00000278267 0.2236183 0.2193012 0.4349965 0.0000000 0.21942369 0.0000000
                   HEA061    HEA062    HEA063     HEA064     HEA065     HEA066
ENSG00000227232 0.3528731 0.0000000 0.0000000 0.06183451 0.10298192 0.05475139
ENSG00000278267 0.2352487 0.2880235 0.1926772 0.49467608 0.15447288 0.49276252
                    HEA067     HEA068     HEA069     HEA070    HEA071
ENSG00000227232 0.00000000 0.06713275 0.06693886 0.04887949 0.1122282
ENSG00000278267 0.42685808 0.13426549 0.46857204 0.09775898 0.3366846
                    HEA181     HEA182     HEA183    HEA184    HEA185     HEA186
ENSG00000227232 0.20035083 0.28009698 0.05465108 0.0000000 0.0000000 0.41517606
ENSG00000278267 0.20035083 0.14004849 0.32790646 0.6985910 0.3036215 0.17793260
 [ reached getOption("max.print") -- omitted 4 rows ]

3.3 Uso de clases específicas para manejar los datos

Cuando se trabaja con distintos objetos referidos a unos mismos datos, como la matriz de contajes y el objeto “targets”, es útil disponer de clases contenedoras que permitan trabajar con todos ellos a la vez, lo que no sólo facilita el trabajo sino que ayuda a evitar “desincronizaciones”.

Éste es el caso de la clase ExpressionSet habitualmente utilizada con microarrays o de la clase que la generaliza, llamada SummarizedExperiment.

Para datos de contaje es habitual usar una clase similar a ExpressionSet llamada DGEList” pensadas para manejar datos de contajes , definida en el paquete edgeR. Esta clase, más simple que las anteriores, utiliza listas para almacenar recuentos de “reads” e información asociada de tecnologías de secuenciación o expresión génica digital. Puede encontrarse información al respecto en la ayuda del paquete edgeR.

An object of class "DGEList"
$counts
                    COV145     COV147    COV149     COV150     COV151
ENSG00000227232 0.06644967 0.00000000 0.2218932 0.84819243 0.07320991
ENSG00000278267 0.19934902 0.07170342 0.1479288 0.00000000 0.21962973
                    COV152    COV153     COV154   COV155    COV156     COV157
ENSG00000227232 0.15060933 1.3880853 0.41957493 1.163313 0.0000000 0.11736232
ENSG00000278267 0.00000000 0.0000000 0.06992916 0.174497 0.1279917 0.11736232
                   COV175    COV176    COV177    COV178     COV179    COV180
ENSG00000227232 0.0000000 0.9503052 1.9885554 0.1726206 0.14628246 0.6558384
ENSG00000278267 0.2236183 0.2193012 0.4349965 0.0000000 0.21942369 0.0000000
                   HEA061    HEA062    HEA063     HEA064     HEA065     HEA066
ENSG00000227232 0.3528731 0.0000000 0.0000000 0.06183451 0.10298192 0.05475139
ENSG00000278267 0.2352487 0.2880235 0.1926772 0.49467608 0.15447288 0.49276252
                   HEA067     HEA068     HEA069     HEA070    HEA071    HEA181
ENSG00000227232 0.0000000 0.06713275 0.06693886 0.04887949 0.1122282 0.2003508
ENSG00000278267 0.4268581 0.13426549 0.46857204 0.09775898 0.3366846 0.2003508
                    HEA182     HEA183    HEA184    HEA185     HEA186
ENSG00000227232 0.28009698 0.05465108 0.0000000 0.0000000 0.41517606
ENSG00000278267 0.14004849 0.32790646 0.6985910 0.3036215 0.17793260
 [ reached getOption("max.print") -- omitted 3 rows ]
33034 more rows ...

$samples
       group lib.size norm.factors samples group.1 cols
COV145 COVID 999825.8            1  COV145   COVID  red
COV147 COVID 999530.1            1  COV147   COVID  red
COV149 COVID 999667.8            1  COV149   COVID  red
COV150 COVID 999872.5            1  COV150   COVID  red
COV151 COVID 999904.2            1  COV151   COVID  red
29 more rows ...

$genes
                          genes
ENSG00000227232 ENSG00000227232
ENSG00000278267 ENSG00000278267
ENSG00000233750 ENSG00000233750
ENSG00000268903 ENSG00000268903
ENSG00000269981 ENSG00000269981
33034 more rows ...

Uno de los aspectos interesantes de estas clases es la posibilidad de extraer partes de todos los objetos a la vez con el operador de “subsetting”.

[1] 33039    34
       group lib.size norm.factors samples group.1 cols
COV145 COVID 999825.8            1  COV145   COVID  red
COV147 COVID 999530.1            1  COV147   COVID  red
COV149 COVID 999667.8            1  COV149   COVID  red
COV150 COVID 999872.5            1  COV150   COVID  red
COV151 COVID 999904.2            1  COV151   COVID  red
COV157 COVID 999582.4            1  COV157   COVID  red
COV175 COVID 998800.7            1  COV175   COVID  red
COV176 COVID 999667.6            1  COV176   COVID  red
COV177 COVID 999807.4            1  COV177   COVID  red
COV178 COVID 999721.9            1  COV178   COVID  red
 [1] "COV145" "COV147" "COV149" "COV150" "COV151" "COV157" "COV175" "COV176"
 [9] "COV177" "COV178"

Aunque podríamos haber creado el objeto a partir de todas las muestras, y haber realizado la extracción de genes y muestras posteriormente, hemos optado por no hacerlo para facilitar el seguimiento del proceso.

3.4 Normalización

Además de estandarizar los contajes, es importante eliminar otros los sesgos de composición entre librerías. Esto puede hacerse aplicando la normalización por el método TMM que genera un conjunto de factores de normalización, donde el producto de estos factores y los tamaños de librería definen el tamaño efectivo de la biblioteca.

La función calcNormFactors, de la librería edgeR, calcula los factores de normalización entre librerías.

Esto no modificará la matriz de contajes, pero actualizará los factores de normalización en el objeto DGEList (sus valores predeterminados son 1).

Es decir, aunque no se observen cambios en la matriz de contajes, cuando se utilizan estos factores de normalización en algún cálculo la importancia de las distintas columnas se tendrá en cuenta.

Resumiendo

Los análisis que se realicen a continuación se basaran en la matriz de contajes, filtrada, estandarizada y normalizada, sobre la que además se toman logaritmo base dos.

Esta será nuestra matriz de partida para los análisis siguientes,

4 Exploración de los datos

Una vez descartados los genes poco expresados y con los recuentos almacenados en un objeto DGEList, podemos`proceder a realizar algunos gráficos exploratorios para determinar si los datos aparentan buena calidad y/o si presentan algun problema.

4.1 Distribución de los contajes

Un diagrama de cajas con los datos, normalizados o no, muestra que la distribución de los contajes es muy asimétrica, lo que justifica la decisión de trabajar con los logaritmos de los datos.

La transformación logarítmica puede hacerse directamente pero es mejor usar la función cpm, como se ha hecho, que agrega una pequeña cantidad para evitar tomar logaritmos de cero.

4.2 Análisis de similaridad entre las muestras

4.2.1 Distancia entre muestras

La función dist permite calcular una matriz de distancias que contiene las comparaciones dos a dos entre todas las muestras. Por defecto se utiliza una distancia euclídea.

       COV145 COV147 COV149 COV150 COV151 COV152 COV153 COV154 COV155 COV156
COV147   80.0                                                               
COV149   70.7   82.0                                                        
       COV157 COV175 COV176 COV177 COV178 COV179 COV180 HEA061 HEA062 HEA063
COV147                                                                      
COV149                                                                      
       HEA064 HEA065 HEA066 HEA067 HEA068 HEA069 HEA070 HEA071 HEA181 HEA182
COV147                                                                      
COV149                                                                      
       HEA183 HEA184 HEA185
COV147                     
COV149                     
 [ reached getOption("max.print") -- omitted 31 rows ]

Las matrices de distancias se pueden visualizar directamente con un heatmap.

Como puede verse _las muestras tienden a agruparse por el factor SANO/COVID, aunque una de las muestras. COV155 se separa del resto de las del grupo COVID.

4.2.2 Agrupamiento jerárquico

Un agrupamiento jerárquico proporciona una representación alternativa, también basada en la matriz de distancias.

Una de las muestras COVID parece más similar a las saludables que a las otras COVID.

4.2.3 Visualización en dimensión reducida

Como puede verse, el gráfico muestra la misma agrupación “natural” y el mismo comportamiento atípico de una muestra.

5 Análisis de expresión diferencial

El objetivo del análisis de expresión diferencial es seleccionar genes cuya expresión difiere entre grupos.

5.1 Selección de genes usando limma-Voom

La ventaja principal de esta aproximación es que permite trabajar con toda la flexibilidad de los modelos lineales para representar diseños experimentales, y, en muchos casos , aprovechar la experiencia previa del usuario en el manejo de limma.

5.1.1 Matriz de diseño y de contrastes

Utilizando la variable group podemos definir una matriz de diseño y, sobre ésta, los contrastes que nos interesan.

       COVID SANO
COV145     1    0
COV147     1    0
COV149     1    0
COV150     1    0
COV151     1    0
COV152     1    0
COV153     1    0
COV154     1    0
COV155     1    0
COV156     1    0
COV157     1    0
COV175     1    0
COV176     1    0
COV177     1    0
COV178     1    0
COV179     1    0
COV180     1    0
HEA061     0    1
HEA062     0    1
HEA063     0    1
HEA064     0    1
HEA065     0    1
HEA066     0    1
HEA067     0    1
HEA068     0    1
HEA069     0    1
HEA070     0    1
HEA071     0    1
HEA181     0    1
HEA182     0    1
HEA183     0    1
HEA184     0    1
HEA185     0    1
HEA186     0    1
attr(,"assign")
[1] 1 1
attr(,"contrasts")
attr(,"contrasts")$group
[1] "contr.treatment"

Dado que estamos interesados en las diferencias entre los grupos, necesitamos especificar qué comparaciones queremos llevar a cabo. Las comparaciones de interés se puede especificar utilizando la función makeContrasts. La matriz de contraste indica qué columnas de la matriz design vamos a comparar. En este caso tan sólo se llevará a cabo una comparación.

       Contrasts
Levels  CONTROLvsCOVID
  COVID              1
  SANO              -1

5.1.2 Transformación de los contajes

An object of class "EList"
$genes
                          genes
ENSG00000227232 ENSG00000227232
ENSG00000278267 ENSG00000278267
ENSG00000233750 ENSG00000233750
ENSG00000268903 ENSG00000268903
ENSG00000269981 ENSG00000269981
33034 more rows ...

$targets
       group  lib.size norm.factors samples group.1 cols
COV145 COVID 1001653.3    1.0018278  COV145   COVID  red
COV147 COVID 1100529.0    1.1010464  COV147   COVID  red
COV149 COVID  955030.2    0.9553476  COV149   COVID  red
COV150 COVID  947796.6    0.9479175  COV150   COVID  red
COV151 COVID  899521.0    0.8996072  COV151   COVID  red
29 more rows ...

$E
                    COV145     COV147     COV149      COV150     COV151
ENSG00000227232 -0.8223650 -1.1381984 -0.4037625  0.50837551 -0.6500950
ENSG00000278267 -0.5183002 -0.9448596 -0.5597126 -0.92265092 -0.3219038
                    COV152     COV153     COV154     COV155     COV156
ENSG00000227232 -0.6246522  0.9358596  0.1073289  1.0159241 -0.8282289
ENSG00000278267 -1.0045156 -0.9810644 -0.5828556 -0.2862517 -0.4994115
                     COV157     COV175     COV176     COV177     COV178
ENSG00000227232 -0.74055345 -0.9859431  0.5497283  1.3062240 -0.6882933
ENSG00000278267 -0.74055345 -0.4526424 -0.4619603 -0.1060517 -1.1161581
                      COV179     COV180      HEA061      HEA062     HEA063
ENSG00000227232 -0.654532686  0.1342634 -0.20656331 -1.01143480 -1.1012204
ENSG00000278267 -0.499855856 -1.0746764 -0.42066201 -0.35512422 -0.6309653
                    HEA064     HEA065      HEA066     HEA067     HEA068
ENSG00000227232 -0.9265423 -0.7259291 -0.78574256 -1.1218171 -0.9590727
ENSG00000278267 -0.1024608 -0.6077105  0.05386471 -0.2313968 -0.7976723
                    HEA069     HEA070     HEA071     HEA181     HEA182
ENSG00000227232 -1.0099326 -0.8215555 -0.6475041 -0.4772445 -0.4479095
ENSG00000278267 -0.2372664 -0.6984811 -0.1968898 -0.4772445 -0.7333818
                    HEA183     HEA184     HEA185     HEA186
ENSG00000227232 -0.8517698 -1.0495769 -0.9932084 -0.1628122
ENSG00000278267 -0.2738825  0.2117625 -0.3086203 -0.5957197
 [ reached getOption("max.print") -- omitted 3 rows ]
33034 more rows ...

$weights
          [,1]      [,2]      [,3]      [,4]      [,5]      [,6]      [,7]
[1,]  6.976324  6.423290  7.345017  7.414329  7.913732  6.966953  7.070951
[2,] 13.344117  9.792171 16.233083 16.785098 21.291716 13.276959 14.091120
          [,8]      [,9]     [,10]     [,11]     [,12]     [,13]     [,14]
[1,]  8.568175  9.129126  8.046130  6.792973  7.049157  7.046103  6.946915
[2,] 27.856177 34.175509 22.694541 12.084889 13.891449 13.863695 13.134382
         [,15]     [,16]     [,17]     [,18]     [,19]     [,20]     [,21]
[1,]  6.497274  6.878669  6.667025 21.376991 19.139568 14.546448 14.827401
[2,] 10.280091 12.659254 11.283734  8.528972  8.192672  7.560903  7.603134
         [,22]     [,23]     [,24]     [,25]     [,26]     [,27]     [,28]
[1,] 20.095913 24.685105 13.715056 13.110198 11.661751 22.943084 24.332783
[2,]  8.340050  8.955089  7.428218  7.308304  7.031769  8.740787  8.912791
         [,29]     [,30]     [,31]     [,32]     [,33]     [,34]
[1,] 22.360747 15.055081 19.758022 16.994284 20.284112 17.780471
[2,]  8.666115  7.636840  8.288627  7.906856  8.368394  8.008534
 [ reached getOption("max.print") -- omitted 3 rows ]
33034 more rows ...

$design
       COVID SANO
COV145     1    0
COV147     1    0
COV149     1    0
COV150     1    0
COV151     1    0
29 more rows ...

5.1.3 Selección de genes diferencialmente expresados

Como en el caso de los microarrays el objeto voomObj y las matrices de diseño y contrastes se utilizaran para ajustar un modelo y, a continuación realizar las comparaciones especificadas sobre el modelo ajustado. El proceso finaliza con la regularización del estimador del error usando la función eBayes.

5.1.4 Top tables

Los resultados de un análisis de expresión diferencial se pueden extraer con la función topTable. Esta función genera una tabla de resultados cuyas columnas contienen información acerca de los genes y la diferencia entre los grupos comparados. Concretamente:

                          genes    logFC  AveExpr        t      P.Value
ENSG00000167900 ENSG00000167900 2.832423 4.111065 12.36669 8.088344e-15
ENSG00000090889 ENSG00000090889 2.476447 1.292066 11.83766 3.015984e-14
ENSG00000111206 ENSG00000111206 2.423923 2.018435 11.80083 3.309678e-14
ENSG00000166851 ENSG00000166851 3.223189 2.615540 11.65570 4.781203e-14
ENSG00000135476 ENSG00000135476 2.236975 1.705254 11.53985 6.425179e-14
ENSG00000123485 ENSG00000123485 2.553959 1.914888 11.46276 7.828950e-14
                   adj.P.Val        B
ENSG00000167900 2.672308e-10 23.54910
ENSG00000090889 3.644948e-10 22.09352
ENSG00000111206 3.644948e-10 21.99541
ENSG00000166851 3.949154e-10 21.70990
ENSG00000135476 4.013294e-10 21.35292
ENSG00000123485 4.013294e-10 21.18953

5.1.5 Visualización de los resultados

Para visualizar los resultados podemos usar un volcanoPlot:

Con el fin de observar si existen perfiles de expresión diferenciados podemo realizar un mapa de colores con los genes más diferencialmente expresados.

Es decir, fijamos un criterio de selección de genes y retenemos aquellos componentes de la tabla de resultados que lo cumplen. Por ejemplo: Genes con un p-balor ajustado inferior a 0.001 y un `fold-change’ superior a 2.

[1] 199

Con la matriz de expresión de los genes que verifican dicha condición se puede construir un heatmap.

Los dos grupos estan diferenciados, sobre todo en un subconjunto de genes en donde las expresiones toman signos (colores) distintos. En otro de los grupos parece que, en el grupo de individuos sanos los genes diferencialmente expresados se encuentran sobre-expresados en el grupo COVID, y apenas expresados en el grupo sanos.

5.2 Análisis de expresión diferencial usando el paquete edgeR

El análisis con edgeR es similar al anterior (se originan en el mismo equipo de investigación) pero la modelización es distinta.

El análisis utiliza un GLM pero, en una forma que recuerda a lo que se hace con limma, realiza un paso adicional,en el que se calcula una estimación mejorada de la dispersión (variabilidad) de las muestras que integra las estimaciones individuales y la global mediante estimación Bayes empírica.

Con este objeto, que añade a los contajes normalizados los estimadores mejorados de dispersión, se ajusta un modelo lineal generalizado con distribución binomial para los errores.

Una vez ajustado el modelo se procede a construir el contraste y realizar el test.

De hecho podemos usar la matriz de contrastes que construímos para limma voom, en la misma forma que hemos reutilizado la de diseño.

An object of class "DGELRT"
$coefficients
                    COVID      SANO
ENSG00000227232 -14.28149 -15.28977
ENSG00000278267 -15.16817 -14.68027
ENSG00000233750 -15.61279 -15.65247
ENSG00000268903 -13.79667 -14.19322
ENSG00000269981 -14.09073 -14.31982
ENSG00000241860 -14.96062 -14.97429

$fitted.values
                    COV145     COV147     COV149     COV150     COV151
ENSG00000227232 0.50362351 0.55333742 0.48018176 0.47654477 0.45227218
ENSG00000278267 0.13404573 0.14727771 0.12780641 0.12683838 0.12037793
                    COV152     COV153     COV154     COV155    COV156
ENSG00000227232 0.50436793 0.49623562 0.42920617 0.41355987 0.4463544
ENSG00000278267 0.13424387 0.13207935 0.11423862 0.11007416 0.1188028
                    COV157     COV175     COV176     COV177    COV178
ENSG00000227232 0.51862932 0.49791657 0.49815307 0.50596775 0.5449482
ENSG00000278267 0.13803971 0.13252676 0.13258970 0.13466968 0.1450448
                    COV179     COV180     HEA061    HEA062     HEA063
ENSG00000227232 0.51149863 0.52950246 0.10257076 0.1050509 0.11179640
ENSG00000278267 0.13614179 0.14093373 0.29173133 0.2987852 0.31797088
                    HEA064     HEA065     HEA066     HEA067     HEA068
ENSG00000227232 0.11129685 0.10394130 0.09967537 0.11340392 0.11490841
ENSG00000278267 0.31655006 0.29562943 0.28349630 0.32254297 0.32682204
                    HEA069     HEA070     HEA071     HEA181     HEA182
ENSG00000227232 0.11899089 0.10109909 0.09995140 0.10161035 0.11090222
ENSG00000278267 0.33843341 0.28754563 0.28428137 0.28899975 0.31542764
                    HEA183     HEA184     HEA185     HEA186
ENSG00000227232 0.10432431 0.10786525 0.10373204 0.10677595
ENSG00000278267 0.29671878 0.30678991 0.29503426 0.30369172
 [ reached getOption("max.print") -- omitted 4 rows ]

$deviance
ENSG00000227232 ENSG00000278267 ENSG00000233750 ENSG00000268903 ENSG00000269981 
      13.203208        3.466706        1.937466        6.511782       10.312153 
ENSG00000241860 
       4.314431 

$method
[1] "oneway"

$unshrunk.coefficients
                    COVID      SANO
ENSG00000227232 -14.50309 -16.07675
ENSG00000278267 -15.82674 -15.03147
ENSG00000233750 -17.00917 -17.18018
ENSG00000268903 -13.92723 -14.39407
ENSG00000269981 -14.27015 -14.55116
ENSG00000241860 -15.45818 -15.48078

$df.residual
[1] 32 32 32 32 32 32

$design
       COVID SANO
COV145     1    0
COV147     1    0
COV149     1    0
COV150     1    0
COV151     1    0
29 more rows ...

$offset
         [,1]    [,2]    [,3]    [,4]     [,5]     [,6]     [,7]     [,8]
[1,] 13.81716 13.9113 13.7695 13.7619 13.70962 13.81864 13.80238 13.65727
         [,9]    [,10]    [,11]    [,12]    [,13]    [,14]    [,15]    [,16]
[1,] 13.62014 13.69645 13.84652 13.80577 13.80624 13.82181 13.89602 13.83268
        [,17]    [,18]    [,19]    [,20]    [,21]    [,22]    [,23]    [,24]
[1,] 13.86727 13.79954 13.82344 13.88567 13.88119 13.81282 13.77091 13.89995
        [,25]    [,26]    [,27]    [,28]    [,29]    [,30]   [,31]    [,32]
[1,] 13.91313 13.94804 13.78509 13.77367 13.79014 13.87764 13.8165 13.84987
       [,33]    [,34]
[1,] 13.8108 13.83972
attr(,"class")
[1] "CompressedMatrix"
attr(,"Dims")
[1]  6 34
attr(,"repeat.row")
[1] TRUE
attr(,"repeat.col")
[1] FALSE

$dispersion
[1] 9.765625e-05 9.765625e-05 9.765625e-05 9.765625e-05 9.765625e-05
[6] 9.765625e-05

$prior.count
[1] 0.125

$AveLogCPM
[1] 1.198866 1.145721 1.023866 1.444349 1.352131 1.128949

$df.residual.zeros
[1] 32 32 32 32 32 32

$df.prior
ENSG00000227232 ENSG00000278267 ENSG00000233750 ENSG00000268903 ENSG00000269981 
       6.412772        9.221763        9.221763        9.221763        9.221763 
ENSG00000241860 
       9.221763 

$var.post
ENSG00000227232 ENSG00000278267 ENSG00000233750 ENSG00000268903 ENSG00000269981 
      0.3701341       0.1170509       0.0646997       0.2054864       0.2930194 
ENSG00000241860 
      0.1356745 

$var.prior
[1] 0.15822624 0.14729701 0.07911392 0.21240299 0.19157096 0.13861893

$samples
       group lib.size norm.factors samples group.1 cols
COV145 COVID 999825.8    1.0018278  COV145   COVID  red
COV147 COVID 999530.1    1.1010464  COV147   COVID  red
COV149 COVID 999667.8    0.9553476  COV149   COVID  red
COV150 COVID 999872.5    0.9479175  COV150   COVID  red
COV151 COVID 999904.2    0.8996072  COV151   COVID  red
29 more rows ...

$genes
                          genes
ENSG00000227232 ENSG00000227232
ENSG00000278267 ENSG00000278267
ENSG00000233750 ENSG00000233750
ENSG00000268903 ENSG00000268903
ENSG00000269981 ENSG00000269981
ENSG00000241860 ENSG00000241860

$table
                      logFC   logCPM            F     PValue
ENSG00000227232  1.45463439 1.198866 13.174669854 0.02722442
ENSG00000278267 -0.70389442 1.145721  9.137260461 0.30104772
ENSG00000233750  0.05724585 1.023866  0.145023951 0.92283279
ENSG00000268903  0.57210237 1.444349  6.405756530 0.25124362
ENSG00000269981  0.33050524 1.352131  1.266644864 0.54237109
ENSG00000241860  0.01973312 1.128949  0.006129263 0.97699496

$comparison
[1] "1*COVID -1*SANO"

$df.test
[1] 1 1 1 1 1 1

$df.total
ENSG00000227232 ENSG00000278267 ENSG00000233750 ENSG00000268903 ENSG00000269981 
       38.41277        41.22176        41.22176        41.22176        41.22176 
ENSG00000241860 
       41.22176 

Los resultados se almacenan en un objeto similar a la ’ topTable’ de limma.

Coefficient:  1*COVID -1*SANO 
                          genes    logFC   logCPM        F       PValue
ENSG00000123485 ENSG00000123485 2.969009 2.874508 143.3978 6.410064e-15
ENSG00000167900 ENSG00000167900 3.052165 4.946915 165.9368 1.024164e-14
ENSG00000187837 ENSG00000187837 1.719557 4.813853 135.1654 1.889383e-14
ENSG00000111206 ENSG00000111206 2.834520 2.919472 132.0747 2.767000e-14
ENSG00000178999 ENSG00000178999 2.987952 3.225967 135.0640 4.056862e-14
ENSG00000166851 ENSG00000166851 3.638997 3.745942 144.7154 5.207210e-14
                         FDR
ENSG00000123485 1.691868e-10
ENSG00000167900 1.691868e-10
ENSG00000187837 2.080777e-10
ENSG00000111206 2.285473e-10
ENSG00000178999 2.680693e-10
ENSG00000166851 2.867350e-10

Podemos seleccionar los genes más diferencialmente expresados de la misma forma que hicimos con limma-voom

[1] 376

Obsérvese que la lista de genes seleccionados es muy similar en ambos casos.

6 Anotación de resultados y análisis de significación biológica

Para el análisis de significación se utilizan dos listas de transcritos:

[1] 383
[1] 33039

6.1 Anotación de los identificadores

Esto es posible, y de hecho sencillo de llevar a cabo, usando el paquete annotate.

 [1] "ACCNUM"       "ALIAS"        "ENSEMBL"      "ENSEMBLPROT"  "ENSEMBLTRANS"
 [6] "ENTREZID"     "ENZYME"       "EVIDENCE"     "EVIDENCEALL"  "GENENAME"    
[11] "GENETYPE"     "GO"           "GOALL"        "IPI"          "MAP"         
[16] "OMIM"         "ONTOLOGY"     "ONTOLOGYALL"  "PATH"         "PFAM"        
[21] "PMID"         "PROSITE"      "REFSEQ"       "SYMBOL"       "UCSCKG"      
[26] "UNIPROT"     
          ENSEMBL SYMBOL ENTREZID                                   GENENAME
1 ENSG00000167900    TK1     7083                         thymidine kinase 1
2 ENSG00000090889  KIF4A    24137                   kinesin family member 4A
3 ENSG00000111206  FOXM1     2305                            forkhead box M1
4 ENSG00000166851   PLK1     5347                         polo like kinase 1
5 ENSG00000135476  ESPL1     9700 extra spindle pole bodies like 1, separase
6 ENSG00000123485  HJURP    55355      Holliday junction recognition protein
[1] 385   4

Como puede verse, el número de anotaciones es el mismo que el de identificadores ENSEMBL, lo que podría llevar a pensar que, es posible que, antes de subir los datos a GEO se hayan agrupado los contajes por genes.

Para la anotación del universo se procederá igual.

          ENSEMBL SYMBOL ENTREZID                                   GENENAME
1 ENSG00000167900    TK1     7083                         thymidine kinase 1
2 ENSG00000090889  KIF4A    24137                   kinesin family member 4A
3 ENSG00000111206  FOXM1     2305                            forkhead box M1
4 ENSG00000166851   PLK1     5347                         polo like kinase 1
5 ENSG00000135476  ESPL1     9700 extra spindle pole bodies like 1, separase
6 ENSG00000123485  HJURP    55355      Holliday junction recognition protein
[1] 33202     4

En este caso se observa como hay más anotaciones que transcritos, lo que sugiere que múltiples transcritos han sido mapeados en el mismo gen.

6.2 Análisis de enriquecimiento

El paquete clusterProfiler admite identificadores de tipo ENSEMBL y permite gran variedad de análisis complementarios al enriquecimiento por lo que, es una de las mejores opciones para el análisis de significación biológica.

                   ID
GO:0006958 GO:0006958
GO:0002455 GO:0002455
GO:0006956 GO:0006956
GO:0016064 GO:0016064
GO:0019724 GO:0019724
                                                              Description
GO:0006958                       complement activation, classical pathway
GO:0002455 humoral immune response mediated by circulating immunoglobulin
GO:0006956                                          complement activation
GO:0016064                        immunoglobulin mediated immune response
GO:0019724                                       B cell mediated immunity
           GeneRatio   BgRatio        pvalue      p.adjust        qvalue Count
GO:0006958    80/342 120/16130 4.046401e-107 1.215134e-103 1.117659e-103    80
GO:0002455    81/342 135/16130 9.256747e-103  1.389901e-99  1.278405e-99    81
GO:0006956    81/342 148/16130  4.506649e-98  4.511156e-95  4.149280e-95    81
GO:0016064    81/342 200/16130  3.900362e-84  2.928197e-81  2.693303e-81    81
GO:0019724    81/342 203/16130  1.728105e-83  1.037900e-80  9.546418e-81    81

Con los resultados del análisis de enriquecimiento se pueden llevar a cabo distintas visualizaciones cuya interpretación exacta puede verse en el manual de clusterProfiler